Vanishing of the negative-sign problem of quantum Monte Carlo simulations in 

one-dimensional frustrated spin systems 
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The negative-sign problem in one-dimensional frustrated quantum spin systems is solved. We 
can remove negative signs of the local Boltzmann weights by using a dimer basis that has the spin- 
reversal symmetry. Validity of this new basis is checked in a general frustrated double-spin-chain 
system, namely the J0-J1-J2-J3 model. The negative sign vanishes perfectly for Jo + Ji < J3. 
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The class of low-dimensional quantum spin systems is 
now enjoying a revival of interest both theoretically and 
experimentally. This recent interest originates in the pos- 
sibility of the superconductivity upon doping a carrier 
to an insulator that has a spin gap above the ground 
state. Spin-ladder models, or more generally double-spin- 
chain models, are known as such candidates. ^ Synthe- 
ses of various corresponding materials [^-^ support the 
progress of this field under the cooperations between the 
experiment and the theory. For example, magnetic sus- 
ceptibility measurements on KCuCls ^ and on CaV205 
Q indicate a spin gap behavior. These systems are con- 
sidered to be explained by a frustrated double-spin-chain 
model, and the strength of each interaction bond can 
be estimated by comparing the experimental data with 
the theoretical ones. [|| In such analyses, we must calcu- 
late the observable quantities for given interaction bonds. 
This is generally a very difficult task. 

Numerical investigations may serve as a powerful tool 
of calculating the thermodynamic quantities, unless the 
problem is solved exactly. For example, we can obtain 
all the eigenvalues of a finite system by the numerical di- 
agonalization. Then any quantity at any temperature is 
calculated, but the size of the system is very restricted; 
up to 16 or 18 of S" = 1/2 spin sites. The quantum Monte 
Carlo (QMC) method can handle much larger systems, if 
there is no frustration. |^ Contrary to this, the sampling 
ruins at low temperatures by the negative-sign problem 
in the frustrated systems. The transition probability may 
take a negative value in this situation. One uses its ab- 
solute value for the update, and then reweights both the 
number of steps and the physical quantity in order to 
obtain the correct value. The problem is caused by the 
fact that the partition function defined by the absolute 
values of the local weights is far from the original one, 
and thus it becomes serious exponentially with the prod- 
uct of the inverse temperature (3 and the system size N. 
We are necessary to overcome this negative-sign problem 
in order to get meaningful numerical data that can be 
compared directly with the experiment. This is the main 
subject of this letter. 

There have been already proposed several techniques 



that relax the sign problem: the auxiliary field method, 
the transfer-matrix Monte Carlo method, j|] the 
reweighting method, ||^ and the restructuring method. 
JlOt However, it remains. Therefore, at low temperatures 
or for large system sizes, the QMC samplings any way 
break down. 

In this letter we demonstrate that the negative-sign 
problem is totally removed for the first time in the non- 
trivial frustrated spin systems. The key idea is an ex- 
tension of the restructuring method, ]lO| combined with 
the spin-reversal symmetry that the system usually pos- 
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FIG. 1. Shape of the general double-spin-chain model we 
treat in this letter. 

We consider the generalized double-spin-chain system 
defined by its next-nearest-neighbor interactions, Jo and 
Ji, and by the alternating nearest-neighbor interactions, 
J2 and J3, as depicted in Fig. |l|. This model can de- 
scribe various systems; it is reduced to the Majumdar- 
Ghosh model |ll| with a choice of the parameter set 
(Jo, Ji, ^2,^3) = (0.5,0.5,1,1). The dimer-fluid tran- 
sition point, where the spin-Peierls material CuGeOa 
]l2[ is suggested to realize, 14| is ( Jq, Ji, J2, J3) = 
(0.2411,0.2411,1,1). |l| 

We divide the above Hamiltonian into two parts for 
the Suzuki- Trotter decomposition. 



(^+1)72 N/2 

Ti-l ^ ^ h2n-l, 'H2 ~ h2n 
n— 1 n—1 

Ki — JqCu ■ Cr-n+l + JlTn ' Tn+1 
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J2 



-^CTn+l ■ T n+1 + Js'T n ' Cn+l, (2) 



Then, the total Boltzmann weight, {4>\ exp[— /STCWtp) , is 
decomposed into the product of the local weights as 
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~\e --\^) = (V'|[e-^«ie-^«=]'"|i/.). 
We propose the representation basis \'ip) in the foUow- 
ing way. Usually, diagonal representation of each spin 
is chosen for Here, we couple two spins, cr„ and 

T„, and consider these two spins as an unit. This dimer 
unit takes four states associated with the eigenvalues 
of each spin, ,t^). We restructure these four states so 
that they become eigenstates of the spin-reversal opera- 
tion R: 



vi = {\ T,T> + 
«2 = (| T,T>- 
^^3 = (| Ta) + 
^'4 = (| T,i>- 



U,i»/V2, 
U,i»/V2, 
U,T))/x/2, 
U,T»/V2. 



(3) 
(4) 
(5) 
(6) 



Here, t and [ denote the eigenstates. Direct prod- 
uct of these four states spans the whole phase space. In 
this new basis, three interaction bonds, Jq, Ji and J3 be- 
come a single effective bond connecting the neighboring 
dimer units, and the J2 bonds only contribute to the in- 
ner energy of the dimer units. It is quite striking that 
this basis transformation alone removes the negative-sign 
problem. Let us call this basis 'dimer-R' basis hereafter 
for simplicity. 

In order to calculate the local plaquette Boltzmann 
weight, {vi,Vj\ exp[— /3/i„/m]|w^,w^), we need a matrix of 
the local Hamiltonian, {vi,Vj\hn\v[,v'j) . Since the Trot- 
ter slice (3/m should be small, the sign of the local Boltz- 
mann weight is known from the sign of the local Hamilto- 
nian as exp[— /3/i„/m] ~ 1 — /3/i„/m. Thus we first write 
down the matrix element of /i„. This 16 x 16 square ma- 
trix is block-diagonalized by the spin-reversal symmetry 
into four 4x4 sub-blocks: 
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where a = J0 + J1 + J3, b ^ -J0 + J1 + J3, c = Jq — Ji + J^, 
and d = — Jo — Ji + J3 . 

These matrix elements include negative signs and it 
seems that the negative-sign problem exists. However, 
we can remove the negative sign by the following nonlo- 
cal unitary transformation. [^^|l8| This transformation 
is an adaptation of the Kennedy- Tasaki transformation 
||l6| of the 5 = 1 AF Heisenberg chain to the S = 1/2 
double-spin-chain systems. It is defined by U : 

N 



where Sn ~ cTn 
transformed as. 
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Then, the local Hamiltonian ft,„ is 
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Each matrix element of this transformed Hamiltonian be- 
comes the one whose negative sign is taken away from 
the original one. Note that the signs of the interactions, 
Jo, Ji, and J3 are changed. All the matrix elements be- 
come positive if all the a, 5, c, and d are positive. This 
is when 



d = - Jo - Ji + J3 > 0, 



(10) 



for the AF positive values of Jq, Ji and J3. This condi- 
tion includes various interesting models as stated before. 
The above condition does not restrict the value of J2, 
since it only contributes to the diagonal matrix elements. 
We can now perform quantum Monte Carlo simulations 
without the negative-sign problem. Before the demon- 
stration of the numerical results, we point out several 
important notices. 

Since we have removed the negative-sign problem by 
using the spin-reversal symmetry, we have to be care- 
ful about the operation that breaks this symmetry. For 
example, the negative-sign problem appears again, if 
we apply the uniform magnetic field. In our new rep- 
resentation, the matrix element for the uniform mag- 
netic field H only takes a value between vi and V2, i.e., 
- HS^\v2) -iJ, for = af + rf . Thus this term 
has a value in the off- block-diagonal part. 
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Another point is that we cannot observe the uniform 
magnetic susceptibihty from the fluctuation of the mag- 
netization, since it does not fluctuate but always remains 
zero. In other words, a local expectation value of the 
magnetization, {vi, [(S'f exp[— /3/ii/m]|w^, v'j), van- 
ishes for a non- vanishing Boltzmann weight, and vice 
versa. Therefore, we have considered the following two 
different methods of calculating the susceptibility. Both 
methods intrinsically suffer the negative-sign problem, 
but, it is not so severe that we can obtain sufficient data 
for the large system sizes at low temperatures. 

The first one is that we take all the configurational 
summation over a single Trotter layer ip' to deduce the 
expectation value of {i{j\{J2i Sf)^ exp[—PHi/m]\i{}') in 
the measurement stage at each Monte Carlo step. That 
is, we take a Monte Carlo summation of 



point. We first check our simulation giving a correct val- 
ues by comparing with the exact results obtained by the 
numerical diagonalization in the system with 10 spins. 
Then we show the results of 34 spins, which cannot be 
achieved by any other method. The boundary conditions 
are set open. The temperatures that we actually ran sim- 
ulations are pointed by arrows in Fig. I and Fig. S The 
others are estimated by the reweighting method. B We 
have done the Trotter extrapolations by using five or six 
different Trotter numbers ranging from P/m — 0.5 to 
about 0.2. The number of the Monte Carlo steps is typ- 
ically five millions divided into ten parts to estimate the 
deviations. The correlation time for the energy is about 
5 steps for T/J2 = 0.2 and m = 24; that for the sus- 
ceptibility obtained by the transfer matrix as mentioned 
above is about 1 step. 
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as a contribution to the numerator of the susceptibility. 
Here, '00, ipi, and tp2 are the states that are actually re- 
alized in the simulation, and ip'^ is the virtual state that 
should be traced out. This tracing-out can be done by 
multiplying a 4 x 4 transfer matrix along the real space 
direction. The contribution to the denominator. 



( V'o I e-'^^i I V'l ) (-01 1 e-''^^ I V'2 
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always becomes unity by the conservation law of the spin- 
reversal eigenvalue and thus it becomes the number of 
steps. This way of calculating the susceptibility gives the 
correct values with sufficiently small error bars especially 
at rather high temperatures. The sampling suddenly be- 
comes worse at low temperatures in larger sizes. This 
is because we simulate the full system and observe the 
susceptibility in the system that a single layer is traced 
out. Since both systems are different from each other 
by one Trotter layer, there should exist the negative-sign 
problem. This appears in both the numerator and the 
denominator of Eq. (p4|), which can take both a positive 
and a negative sign. 

The other method is the numerical differentiation of 
the magnetization when we apply the sufficiently small 
magnetic field. As discussed above, the negative-sign 
problem appears under the magnetic field, however, the 
problem is not so serious if it is small enough. For ex- 
ample, the typical negative-sign ratio at the temperature 
T/J2 = 0.04 is only 0.5 with the field H/J2 = 0.02 in 
the system with 34 spins. Thus we can calculate the 
susceptibility without any difficulty. 

Now that we show the numerical results and demon- 
strate the usefulness of the present new basis. We set 
the interaction bonds Jq = Ji and J2 = J3 with J1/J2 = 
0.2411, so that this system is at the dimer- fluid transition 
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FIG. 2. (a) The susceptibility (circles) and the specific heat 
(triangles) for the system with 10 spins calculated by using 
the present dimer-R basis. The exact results obtained by the 
numerical diagonalization are denoted by lines, (b) Those for 
the same system as (a) but calculated by using the conven- 
tional basis. Arrows indicate the temperatures that the 
simulations are actually performed. The others are estimated 
by the reweighting method. 



Figure ^ (a) shows the susceptibility and the specific 
heat for the system with 10 spins calculated by using 
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the present dimer-R basis; (b) shows the results for the 
same systera calculated by the same program but only the 
representation basis is different using the conventional 
one. The exact results obtained by the numerical diago- 
nalization are plotted by lines. Agreements of the data in 
Fig. H (a) with the exact values are very excellent down 
to the temperature T/ J2 — 0.02 for both quantities. Er- 
ror bars are mostly as small as invisible. We calculated 
the susceptibility here by using the transfer matrix de- 
scribed above as the first method. Since the system size 
is small enough, the instability of the data due to the 
tracing-out of one layer scarcely occurred. On the other 
hand, the samplings of the basis crashes already at 
T/J2 = 0.1. This means that the QMC method cannot 
be applied even in a small system that can be fully diag- 
onalized numerically, once one uses the conventional 
basis. 
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FIG. 3. The susceptibility (circles) and the specific heat 
(triangles) for the system with 34 spins calculated by using 
the present dimer-R basis. Arrows indicate the temperatures 
that the simulations are actually performed. 

Figure H shows the results for the system with 34 spins. 
The susceptibility was calculated by the transfer matrix 
for the data of T/J2 > 0.2, and by the numerical differ- 
entiation of the magnetization for those of T/ J2 < 0.2. 
Both methods give consistent results within the error 
bars, and we have selected the one whose error bar is 
smaller than the other one. The magnitude of the uni- 
form magnetic field for the latter case is H/J2 = 0.02. 
The susceptibility grows as decreasing the temperature, 
i.e., the pseudo gap caused by the finite size of the system 
seems to be much smaller than the temperature range we 
simulated. The specific heat was calculated in the sim- 
ulation with zero field. The error bars for this value is 
also negligible. Rather large ones, though they are within 
the symbols, are solely from the reweighting error. Thus, 
there is no difficulty in calculating the response functions 
of the operators that conserve the spin-reversal symme- 



try. We consider that the temperature can be lowered 
from T/J2 = 0.04, if one can have much more CPU time 
to handle the simulations with large Trotter number. 

We have presented a new possibility of the quantum 
Monte Carlo method by using the representation basis 
that has the spin-reversal symmetry. The negative-sign 
problem is totally removed in the wide range of the pa- 
rameter space of the one-dimensional frustrated spin sys- 
tems. This made possible that one can perform meaning- 
ful simulations within the very restricted computational 
facilities; all the numerical results presented in this letter 
were obtained by a DEC Alpha-433 PC in two weeks. 
These data will not be obtained even by any supercom- 
puter, if one uses the conventional basis. The ex- 
tension of the condition ( p^ ) that the sign problem is 
removed, and also the application to the higher dimen- 
sions are the future problem. In such cases, one should 
take into account the relevant symmetry that the model 
possesses as was successful in the present model. 
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